Assignment 10: Surface currents from the OSCAR model¶
Let's access and plot global surface currents from OSCAR: https://www.esr.org/research/oscar/oscar-surface-currents/ according to Earth and Space Research, who create this data product:
"The Ocean Surface Current Analyses Real-time (OSCAR) is a NASA funded research project and is a direct computation of global surface currents using satellite sea surface height, wind, and temperature. Currents are calculated using a quasi-steady geostrophic model together with an eddy viscosity based wind-driven ageostrophic component and a thermal wind adjustment. The model calculates a surface current averaged over the top 30 m of the upper ocean.
In [1]:
# Our usual imports
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import cool_maps.plot as cplt
import cartopy.crs as ccrs
import cmocean.cm as cmo
import numpy as np
import cartopy.feature as cfeature
%matplotlib inline
The data is available from the NASA Physical Oceanography Distributed Active Archive Center: https://podaac.jpl.nasa.gov/
We can find more information about this data product at https://podaac.jpl.nasa.gov/dataset/OSCAR_L4_OC_third-deg?ids=&values=&search=oscar&provider=PODAAC
Download the 2018 dataset off Canvas here. It used to be accessible at https://podaac-opendap.jpl.nasa.gov/opendap/allData/oscar/preview/L4/oscar_third_deg/oscar_vel2018.nc.gz
In [2]:
# Download this file off Canvas
local = '/Users/csmit/Downloads/oscar_vel2018.nc.gz'
ds = xr.open_dataset(local)
ds
Out[2]:
<xarray.Dataset> Size: 1GB
Dimensions: (time: 72, year: 72, depth: 1, latitude: 481, longitude: 1201)
Coordinates:
* time (time) datetime64[ns] 576B 2018-01-01 2018-01-06 ... 2018-12-26
* year (year) float32 288B 2.018e+03 2.018e+03 ... 2.019e+03 2.019e+03
* depth (depth) float32 4B 15.0
* latitude (latitude) float64 4kB 80.0 79.67 79.33 ... -79.33 -79.67 -80.0
* longitude (longitude) float64 10kB 20.0 20.33 20.67 ... 419.3 419.7 420.0
Data variables:
u (time, depth, latitude, longitude) float64 333MB ...
v (time, depth, latitude, longitude) float64 333MB ...
um (time, depth, latitude, longitude) float64 333MB ...
vm (time, depth, latitude, longitude) float64 333MB ...
Attributes: (12/15)
VARIABLE: Ocean Surface Currents
DATATYPE: 1/72 YEAR Interval
DATASUBTYPE: unfiltered
GEORANGE: 20 to 420 -80 to 80
PERIOD: Jan.01,2018 to Dec.26,2018
year: 2018
... ...
source: Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
contact: Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
company: Earth & Space Research, Seattle, WA
reference: Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
note1: Maximum Mask velocity is the geostrophic component at all...
note2: Longitude extends from 20 E to 420 E to avoid a break in ...Let's first plot global surface U (zonal speed) on a particular time record: 2018-01-01¶
In [3]:
ds.u.sel(time='2018-01-01', method='nearest').plot();
So we can see it is a global dataset. How about we just look at the western north atlantic?¶
Select only lon = 275:350, lat = 60:20, for a single time record: 2018-01-01
remember to use .sel( dim = slice(start, stop)
In [4]:
ds.u.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) ).plot( vmin=-.75 , vmax=.75, cmap='RdBu');
In [5]:
ds.v.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) ).plot( vmin=-.75 , vmax=.75, cmap='RdBu');
Let's save that exact subset you just created into a new dataarray called subset. This will streamline the coding as we move forward, rather than calling the entire dataset each time.¶
In [6]:
subset = ds.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) )
subset
Out[6]:
<xarray.Dataset> Size: 878kB
Dimensions: (year: 72, depth: 1, latitude: 121, longitude: 226)
Coordinates:
time datetime64[ns] 8B 2018-01-01
* year (year) float32 288B 2.018e+03 2.018e+03 ... 2.019e+03 2.019e+03
* depth (depth) float32 4B 15.0
* latitude (latitude) float64 968B 60.0 59.67 59.33 ... 20.67 20.33 20.0
* longitude (longitude) float64 2kB 275.0 275.3 275.7 ... 349.3 349.7 350.0
Data variables:
u (depth, latitude, longitude) float64 219kB ...
v (depth, latitude, longitude) float64 219kB ...
um (depth, latitude, longitude) float64 219kB ...
vm (depth, latitude, longitude) float64 219kB ...
Attributes: (12/15)
VARIABLE: Ocean Surface Currents
DATATYPE: 1/72 YEAR Interval
DATASUBTYPE: unfiltered
GEORANGE: 20 to 420 -80 to 80
PERIOD: Jan.01,2018 to Dec.26,2018
year: 2018
... ...
source: Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
contact: Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
company: Earth & Space Research, Seattle, WA
reference: Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
note1: Maximum Mask velocity is the geostrophic component at all...
note2: Longitude extends from 20 E to 420 E to avoid a break in ...So we can visualize u and v, but how can we combine them? By plotting the vectors!¶
We can use the matplotlib function called .quiver() to plot vectors. It takes as arguments the x- and y- locations, and the u and v components of the velocity: plt.quiver(x, y, u, v)
The velocity components in our data have 3 dims: depth, latitude, and longitude, and quiver can only handle 2D data, so we need to select a dimension using the numpy style selections. In our case the u and v data we want for plotting can be gotten with subset.u[0,:,:] and subset.v[0,:,:]. This is useful notation if the third dimension is more than 1 count, since you can specify which index you want directly without having to re-subset the data.
In [7]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:] );
How about we zoom in, now that we've plotted the full (albeit subsetted) map, using plt.xlim() and plt.ylim()¶
In [8]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:] );
plt.xlim([280, 300]);
plt.ylim([25, 45]);
.quiver works just like .scatter in that you can specify an array for the colors as well¶
To do this, you need to get the magnitude of the velocity. If you picture a single u and v pair as sides of a right triangle: the u is your adjacent side, the v is your opposite side, how do we get the magnitude or the hypotenuse?¶
In [9]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:], (subset.u[0, :, :]**2 + subset.v[0, :, :]**2 )**0.5 )
plt.xlim([280, 300]);
plt.ylim([25, 45]);
plt.colorbar();
Now we know how our data looks, without any bells & whistles. This should always be your first step when working with a new dataset. Plot the data as simply as possible to make sure it looks right, before starting on more complicated visualizations.¶
We've done our data check, so we can move on and plot our data on top of a nice map from Mike Smith's cool_maps. Note cool_maps may struggle with the global visuals, since it's designed for more regional focus.¶
We have to do a bit more work with the latitude and longitude values to make this work.¶
In [10]:
print(subset.longitude.shape,subset.latitude.shape,subset.u.shape,subset.v.shape)
(226,) (121,) (1, 121, 226) (1, 121, 226)
In [11]:
subset = subset.squeeze()
print(subset.longitude.shape,subset.latitude.shape,subset.u.shape,subset.v.shape)
(226,) (121,) (121, 226) (121, 226)
While having that mismatch of shapes works fine for a regular quiver plot because Python can figure out what needs to be matched to make it work, when we plot with cartopy or Mike Smith's cool_maps, we have to use transform. Doing so makes it too hard for Python to guess what needs to be done with the latitude and longitude to make them match the u and v. So we have to be a bit more detailed in our coding.¶
In [12]:
# First we turn lon and lat into variables with shape (121,226).
# The longitude defaulted to 0:360, but cool_maps works with -180:180, so we subtract it by 360.
lon, lat = np.meshgrid(subset.longitude.values - 360,subset.latitude.values)
print('Before: ',lon.shape,lat.shape,subset.u.shape,subset.v.shape)
# Then we have to turn them all into 1D arrays for the transform to work properly. reshape(-1) does just that!
lon = lon.reshape(-1)
lat = lat.reshape(-1)
u = subset.u.values.reshape(-1)
v = subset.v.values.reshape(-1)
umag= (u**2 + v**2 )**0.5
print('After: ',lon.shape,lat.shape,u.shape,v.shape)
Before: (121, 226) (121, 226) (121, 226) (121, 226) After: (27346,) (27346,) (27346,) (27346,)
In [13]:
# Now we set the cool_maps extent
extent = [-80, -40, 20, 45]
# Then borrow our cool_maps create line from past coding exercises
fig,ax = cplt.create(extent, proj=ccrs.Mercator(), figsize=(10,6), oceancolor='w')
ax_posi = plt.quiver(lon,lat,u,v,umag,transform=ccrs.PlateCarree())
# Add in a title
plt.title('OSCAR Surface Currents on January 1, 2018');
# And a colorbar
plt.colorbar(label='Speed [m/s]');
Can we make a streamplot as well?¶
https://matplotlib.org/stable/plot_types/arrays/streamplot.html
We need to np.flip() our latitude, because it goes from 60 N to 20 N, but .streamplot needs values for x and y to be monotonically increasing. Since we np.flip() our latitude, we need to np.flipud() our u and v as well.
In [14]:
lon, lat = np.meshgrid(subset.longitude - 360,np.flip(subset.latitude))
u = np.flipud(subset.u.values.squeeze())
v = np.flipud(subset.v.values.squeeze())
umag= (u**2 + v**2 )**0.5
plt.streamplot(lon,lat,u,v,color=umag,linewidth=umag,density=10,arrowsize=.5);
plt.colorbar();
hh
#2. Plot the annual mean of the surface currents (both global and regional subset) for 2018 as colored vectors.¶
Be sure to check the documentation on xarray means to ensure you take it along the correct dimension and skip NaNs.¶
additional imports¶
In [15]:
# Select surface level (assuming depth=15.0 is the only level)
surface_currents = ds.sel(depth=15.0)
# Calculate the annual mean for 2018
u_mean = surface_currents['u'].mean(dim='time', skipna=True)
v_mean = surface_currents['v'].mean(dim='time', skipna=True)
# Calculate magnitude of the currents
umag_g =(u_mean**2 + v_mean**2)**0.5
# Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines() # Add coastlines
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white') # White background for land
ax.add_feature(cfeature.OCEAN, color='white') # White background for oceans
quiver = plt.quiver(surface_currents['longitude'], surface_currents['latitude'], u_mean, v_mean, umag_g, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='viridis')
ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
plt.title('Annual Mean Surface Currents (Global) - 2018')
plt.show()
# Regional subset for the northwestern Atlantic
region = surface_currents.sel(latitude=slice(45, 25), longitude=slice(280, 300))
u_mean_region = region['u'].mean(dim='time', skipna=True)
v_mean_region = region['v'].mean(dim='time', skipna=True)
umag_r = (u_mean_region**2 + v_mean_region**2)**0.5
# Regional plot
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')
ax.set_xlabel('longitude')
quiver = plt.quiver(region['longitude'], region['latitude'], u_mean_region, v_mean_region, umag_r, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='viridis')
cbar = plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
# Add gridlines with proper longitude and latitude labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
plt.title('Annual Mean Surface Currents (Northwestern Atlantic) - 2018', fontsize=14)
# Adjust layout to prevent squishing
#plt.tight_layout()
plt.show()
#3. Plot the annual mean of the surface currents (both global and regional subset) for 2017 (available off Canvas as well) as colored vectors.¶
In [16]:
local = '/Users/csmit/Downloads/oscar_vel2017.nc.gz'
ds2 = xr.open_dataset(local)
ds2
Out[16]:
<xarray.Dataset> Size: 1GB
Dimensions: (time: 72, year: 72, depth: 1, latitude: 481, longitude: 1201)
Coordinates:
* time (time) datetime64[ns] 576B 2017-01-01 2017-01-06 ... 2017-12-26
* year (year) float32 288B 2.017e+03 2.017e+03 ... 2.018e+03 2.018e+03
* depth (depth) float32 4B 15.0
* latitude (latitude) float64 4kB 80.0 79.67 79.33 ... -79.33 -79.67 -80.0
* longitude (longitude) float64 10kB 20.0 20.33 20.67 ... 419.3 419.7 420.0
Data variables:
u (time, depth, latitude, longitude) float64 333MB ...
v (time, depth, latitude, longitude) float64 333MB ...
um (time, depth, latitude, longitude) float64 333MB ...
vm (time, depth, latitude, longitude) float64 333MB ...
Attributes: (12/15)
VARIABLE: Ocean Surface Currents
DATATYPE: 1/72 YEAR Interval
DATASUBTYPE: unfiltered
GEORANGE: 20 to 420 -80 to 80
PERIOD: Jan.01,2017 to Dec.26,2017
year: 2017
... ...
source: Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
contact: Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
company: Earth & Space Research, Seattle, WA
reference: Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
note1: Maximum Mask velocity is the geostrophic component at all...
note2: Longitude extends from 20 E to 420 E to avoid a break in ...In [18]:
# Select surface level (assuming depth=15.0 is the only level)
surface_currents2 = ds2.sel(depth=15.0)
# Calculate the annual mean for 2018
u_mean2 = surface_currents2['u'].mean(dim='time', skipna=True)
v_mean2 = surface_currents2['v'].mean(dim='time', skipna=True)
# Calculate magnitude of the currents
umag_g2 =(u_mean2**2 + v_mean2**2)**0.5
# Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines() # Add coastlines
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white') # White background for land
ax.add_feature(cfeature.OCEAN, color='white') # White background for oceans
quiver = plt.quiver(surface_currents2['longitude'], surface_currents2['latitude'], u_mean2, v_mean2, umag_g2, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='viridis')
ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
plt.title('Annual Mean Surface Currents (Global) - 2017')
plt.show()
# Regional subset for the northwestern Atlantic
region2 = surface_currents2.sel(latitude=slice(45, 25), longitude=slice(280, 300))
u_mean_region2 = region2['u'].mean(dim='time', skipna=True)
v_mean_region2 = region2['v'].mean(dim='time', skipna=True)
umag_r2 = (u_mean_region**2 + v_mean_region**2)**0.5
# Regional plot
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')
quiver = plt.quiver(region2['longitude'], region2['latitude'], u_mean_region2, v_mean_region2, umag_r2, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='viridis')
cbar = plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
# Add gridlines with proper longitude and latitude labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
plt.title('Annual Mean Surface Currents (Northwestern Atlantic) - 2017')
plt.show()
Noticable change in the regional¶
#4. Plot the difference in annual means of the surface currents (both global and regional subset) between 2017 and 2018 as colored vectors. Describe the differences.¶
You can use plt.clim(vmin,vmax) to set your vmin and vmax outside of the quiver line, to make it easier to center your colorbar on 0. Be sure to pick a colormap that works well for this.¶
In [26]:
surface_currents_2017 = ds2.sel(depth=15.0)
surface_currents_2018 = ds.sel(depth=15.0)
u_mean_2017 = surface_currents_2017['u'].mean(dim='time', skipna=True)
v_mean_2017 = surface_currents_2017['v'].mean(dim='time', skipna=True)
u_diff = u_mean_2018 - u_mean_2017
v_diff = v_mean_2018 - v_mean_2017
# Compute the magnitude of the difference vectors
umag_diff = (u_diff**2 + v_diff**2)**0.5
max_diff = max(np.max(np.abs(u_diff)), np.max(np.abs(v_diff)))
# 1. Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')
# Quiver plot for the global difference
quiver = plt.quiver(surface_currents_2018['longitude'], surface_currents_2018['latitude'], u_diff, v_diff, umag_diff, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='coolwarm', clim=(-max_diff, max_diff))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.colorbar(quiver, label='Surface Current Difference Magnitude (m/s)', orientation='vertical')
# Gridlines with labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# Title
plt.title('Difference in Annual Mean Surface Currents (Global) - 2018 vs 2017', fontsize=14)
plt.show()
# 2. Regional plot (Northwestern Atlantic)
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')
# Quiver plot for the regional difference
quiver = plt.quiver(region_2018['longitude'], region_2018['latitude'], u_diff_region, v_diff_region, umag_diff_r, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='coolwarm', clim=(-max_diff, max_diff))
#Color bar
plt.colorbar(quiver, label='Surface Current Difference Magnitude (m/s)', orientation='vertical')
# Gridlines with labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# Title
plt.title('Difference in Annual Mean Surface Currents (Northwestern Atlantic) - 2018 vs 2017', fontsize=14)
plt.show()
The above plots show the difference between 2018 currents and 2017 currents. The red is positive indicating stronger currents in 2018 comapered to 2017. In most cases curents stayed roughly the same. There is indication that the more prominent currents depicted in both years show a stronger average current in 2018. In the regional model there are red current vectors in in the top right of the plot showing a stronger gulf stream current during 2018.
#5. Plot the annual mean of the surface currents (both global and regional subset) for 2018 as colored streamlines.¶
In [35]:
subset_global = surface_currents_2018
# Meshgrid for longitude and latitude
lon_global, lat_global = np.meshgrid(subset_global.longitude - 360, np.flip(subset_global.latitude))
# Flips u, v for plotting
u_global = np.flipud(u_mean_2018.values)
v_global = np.flipud(v_mean_2018.values)
umag_global = np.sqrt(u_global**2 + v_global**2)
# Global streamline plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
# Coastlines and features
ax.coastlines(resolution='10m', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='lightblue')
# Plot streamlines
strm = plt.streamplot(lon_global, lat_global, u_global, v_global, color=umag_global, linewidth=1, cmap='viridis', density=8, arrowsize=1, transform=ccrs.PlateCarree())
# Color bar
plt.colorbar(strm.lines, label='Surface Current Magnitude (m/s)')
ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
# Title
plt.title('Global Annual Mean Surface Currents (2018)', fontsize=14)
plt.show()
# Regional streamline plot for Northwestern Atlantic
subset_regional = surface_currents_2018.sel(latitude=slice(45, 25), longitude=slice(280, 300))
lon_regional, lat_regional = np.meshgrid(subset_regional.longitude - 360, np.flip(subset_regional.latitude))
# Flip u, v for plotting
u_regional = np.flipud(subset_regional['u'].mean(dim='time', skipna=True).values)
v_regional = np.flipud(subset_regional['v'].mean(dim='time', skipna=True).values)
umag_regional = np.sqrt(u_regional**2 + v_regional**2)
# Regional streamline plot
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280 - 360, 300 - 360, 25, 45], crs=ccrs.PlateCarree())
# Coastlines and features
ax.coastlines(resolution='10m', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='lightblue')
# Plot streamlines
strm = plt.streamplot(lon_regional, lat_regional, u_regional, v_regional, color=umag_regional, linewidth=1, cmap='viridis', density=8, arrowsize=1, transform=ccrs.PlateCarree())
# Color bar
plt.colorbar(strm.lines, label='Surface Current Magnitude (m/s)')
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# Title
plt.title('Regional Annual Mean Surface Currents (Northwestern Atlantic, 2018)', fontsize=14)
plt.show()
In [ ]: